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PREFACE 


The  developments  reported  herein  were  sponsored  by  the  Office, 

Chief  of  Engineers  (OCE),  US  Army,  as  part  of  the  Environmental  and 
Water  Quality  Operational  Studies  (EWQOS)  Program,  under  the  work  unit 
(CWIS  31604)  entitled  "Techniques  to  Meet  Environmental  Quality  Objec¬ 
tives  for  Reservoir  Releases."  The  effort  was  supported  by  Task  IIIA.4, 
"Selective  Withdrawal  Characteristics  of  Various  Outlet  Configurations." 
The  OCE  Technical  Monitors  for  EWQOS  were  Mr.  Earl  E.  Elfcer,  Dr.  John 
Bushman,  and  Mr.  James  L.  Gottesman. 

This  project  was  conducted  during  the  period  June  1979  to  September 
1983  as  part  of  a  contract  (DACW39-78-C-0054)  with  Mississippi  State 
University  (MSU)  to  develop  a  boundary-fitted-coordinate  numerical 
hydrodynamic  code  for  reservoir  selective  withdrawal.  The  research  and 
development  were  carried  out  by  Dr.  Joe  F.  Thompson  of  the  Department  of 
Aerospace  Engineering,  MSU,  and  by  Dr.  Robert  S.  Bernard  of  the  Reser¬ 
voir  Water  Quality  Branch  (RWQB) ,  Hydraulics  Laboratory  (HL) ,  US  Army 
Engineer  Waterways  Experiment  Station  (WES).  Drs.  Thompson  and  Bernard 
prepared  this  report. 

Mr.  H.  B.  Simmons,  Chief,  HL,  and  Mr.  J.  L.  Grace,  Jr.,  Chief  of 
the  Hydraulic  Structures  Division,  directed  the  effort.  Supervision  was 
provided  by  Mr.  J.  P.  Holland,  Chief,  RWQB,  and  Dr.  D.  R.  Smith,  former 
Chief  of  the  RWQB.  Mr.  Mark  S.  Dortch,  formerly  of  the  RWQB,  and 
Dr.  Billy  H.  Johnson  of  the  Mathematical  Modeling  Group,  HL,  monitored 
the  contract.  Program  Manager  of  EWQOS  was  Dr.  J.  L.  Mahloch,  WES. 

During  the  preparation  of  this  report,  COL  Tilford  C.  Creel,  CE, 
and  COL  Robert  C.  Lee,  CE,  were  Commanders  and  Directors  of  WES  and 
Mr.  F.  R.  Brown  was  Technical  Director.  At  the  time  of  publication, 

COL  Allen  F.  Grum,  USA,  was  Director  and  Dr.  Robert  W.  Whalin  waf 
Technical  Director. 


This  report  should  be  cited  as  follows: 


Thompson,  J.  F.,  and  Bernard,  R.  S.  1985.  "WESSEL*. 

Code  for  Numerical  Simulation  of  Two-Dimensional  Time- 
Dependent  Width-Averaged  Flows  with  Arbitrary  Boundaries 
Technical  Report  E-85-8,  US  Army  Engineer  Waterways 
Experiment  Station,  Vicksburg,  Miss. 
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WESSEL:  CODE  FOR  NUMERICAL  SIMULATION  OF  TWO-DIMENSIONAL  TIME- 
DEPENDENT  WIDTH-AVERAGED  FLOWS  WITH  ARBITRARY  BOUNDARIES 


PART  I:  INTRODUCTION 

1.  The  code  WESSEL  solves  the  two  momentum  equations,  the  energy 
equation,  and  the  continuity  equation  on  a  two-dimensional  field  with 
boundaries  of  arbitrary  shape,  including  multiple  inlets,  outlets,  and 
obstacles.  The  basis  of  this  numerical  solution  is  a  boundary-fitted 
curvilinear  coordinate  system  that  allows  all  computation  to  be  done  on 
a  rectangular  field  with  a  square  grid  in  the  transformed  plane,  regard¬ 
less  of  the  boundary  shape  and  configuration  in  the  physical  plane. 

2.  The  finite  difference  solution  is  done  in  finite  volume  formula' 
tion,  as  discussed  in  detail  in  Thompson  and  Bernard  (1985).  The  solu¬ 
tion  is  implicit  in  time,  with  all  the  difference  equations  being  solved 
simultaneously  by  SOR  iteration  at  each  time  step.  The  code  reads  the 
boundary-fitted  coordinate  system  from  the  output  of  the  coordinate  code 
WESCOR,  described  In  Thompson  (1983).  The  input  allows  any  portions  of 
the  boundary  (external  or  obstacles)  to  be  designated  as  inlets,  out¬ 
lets,  no-slip  surfaces,  or  slip  surfaces.  Arbitrary  specification  of 
the  variables  on  inlets  and  outlets  is  allowed.  The  output  is  in  the 
form  of  field  arrays  and  plots  (Figure  1)  of  the  velocity  components, 
pressure,  and  temperature.  All  computation  is  done  in  metric  units,  but 
the  input  and  output  units  may  be  specified  otherwise.  Appendixes  A,  B, 
and  C  contain  input  instructions  for  the  hydrodynamic  code  (WESSEL),  the 
contour  plot  code  (CONTUR) ,  and  the  vector  plot  code  (VECTOR),  respec¬ 
tively.  Sample  runstreams  (input)  are  given  in  Appendix  D. 


PART  II:  CONFIGURATIONS 


Computational  Field 

3.  The  computational  field,  i.e.,  the  transformed  region,  is  com¬ 
posed  of  contiguous  rectangular  blocks  with  a  uniform  square  grid  over 
the  entire  field,  where  the  independent  variables  are  the  curvilinear 
coordinates  (C,n)*  designated  in  the  code  as  the  integers  (I,J).  As 
noted  in  Thompson  and  Bernard  (1985),  the  increments,  AC  and  An  ,  are 
irrelevant  since  they  cancel  out  of  the  difference  equations,  and  hence 
are  made  unity  for  convenience.  The  extent  of  C  and  n  on  the  compu¬ 
tational  field  is  from  1  to  IMAX  and  JMAX,  respectively.  The  dimensions 
of  the  field  arrays  are  (0: IDIM.O: JDIM) .  All  field  arrays  are  extended 
one  line  beyond  the  computational  field  for  convenience.  Thus,  the  max¬ 
imum  values  of  IMAX  and  JMAX  allowed  are  IDIM-1  and  JDIM-1 ,  respectively 

4.  The  coordinate  system  is  generated  with  twice  as  many  points  in 
each  direction  as  are  used  for  the  flow  solution.  Therefore,  the  coordi¬ 
nate  arrays  are  dimensioned  (0: ICIM,0: JCIM) ,  where  ICIM=2*IDIM-2  and 
JClM=2*JDIM-2.  The  point  in  the  coordinate  arrays  corresponding  to  the 
Held  puint  (I,J)  is  thus  (II, JJ),  where  II«2*1-1  and  JJ=2*J-1.  This 
correspondence  is  indicated  in  Figure  2. 


0:  solution  points 

Figure  2.  Diagram  of  the  general  field 


*  For  convenience,  symbols  are  listed  and  defined  in  the  Notation 
(Appendix  K) . 


Point  Classification 


5.  This  code  is  applicable  to  a  wide  variety  of  configurations 
including  multiple  inlets,  outlets,  and  obstacles,  with  all  boundaries 
being  of  arbitrary  shape.  In  order  to  achieve  this  versatility,  it  is 
necessary  that  the  points  be  classified  as  to  location  on  the  various 
types  of  boundaries,  in  the  free  field  or  out  of  the  computation  region. 
This  classification  is  done  by  associating  a  different  name  with  each 
type  of  point  as  follows: 

6.  Points  not  on  boundaries. 

FIELD:  field  point  (not  on  a  boundary) 

•  □ 


0UT:  point  totally  ^ut  of  the  computation  region,  i.e.,  inside 
an  obstacle  or  beyond  the  outer  boundary  of  the  flow 
field 


or 


7 .  Boundary  points  not  on  corners. 

WALB0T:  point  on  a  lower  field  boundary  with  no-slip 


or 


WALT0P:  point  on  a  top  field  boundary  with  no-slip 


WALLEF:  point  on  a  left  field  boundary  with  no-slip 


WALRIT:  point  on  a  right  field  boundary  with  no-slip 


With  the  "WAL"  in  the  above  replaced  by  "SLI,"  the  location  is  the  same, 
but  the  boundary  has  slip.  Similarly  "IN"  in  place  of  "WAL"  refers  to 
inlet,  "0UT"  to  outlet,  and  "SUR"  to  free  surface. 

8.  Corner  points. 

CAVWBL:  point  on  bottom-left  concave  corner  of  field  with 
no-slip 


CAVWTR:  point  on  top-right  concave  corner  of  field  with 
no-slip 


CAVWBR:  point  on  bottom-right  concave  corner  of  field  with 
no-slip 


VEXWBL:  point  on  bottom-left  convex  corner  of  field  with 
no-slip 


VEXWTL:  point  on  top-left  convex  corner  of  field  with  no-slip 


VEXWTR:  point  on  top-right  convex  corner  of  field  with  no-slip 


VEXWBR:  point  on  bottom-right  convex  corner  of  field  with 
no-s lip 


With  the  "W"  in  the  preceding  eight  types  replaced  by  "S,"  the  reference 
is  to  a  slip  boundary. 


9.  Point  notation. 

Each  of  the  point 

type  names  is 

made  an 

integer  and  is  assigned  a 

number.  Thus  0UT-O,  FIELD-25, 

etc.  The  types 

are  also  grouped  as  follows: 

WALB0T=1 

CAVWBL-5 

WALLEF-2 

CAVWTL-6 

WALT0P-3 

CAVWTR-7 

WALRIT-4 

WALL=4 

CAVWBR-8 

WALCAV-8 

VEXWBL=9 

SLIB0T-13 

VEXWTL=10 

SLILEF-14 

VEXWTR- 1 1 

SLIT0P-15 

VEXWBR= 1 2 

WALVEX= 1 2 

SL1RIT-16 

SLIP- 16 

CAVSBL-17 

VEXSBL-21 

CAVSTL-18 

VEXSTL-22 

CAVSTR-19 

VEXSTR-23 

CAVSBR-20 

SLICAV=20 

VEXSBR-24 

SLIVEX-24 

INB0T=26 

0UTB0T-3O 

INLEF-27 

0UTLEF-31 

INT0P=28 

0UTT0P-32 

INRIT=29 

INLET-29 

0UTRIT-33 

0UTLET-33 

SURB0T-34 

SURLEF-35 

SURT0P-36 

SURRIT-37 

SURF-37 

Then  if  a  point  type  is  >0 UT  and  ^WALL,  the  point  must  be  on  one  of  the 
no-slip  walls.  Similarly,  If  a  point  type  is  >  FIELD  and  ^INLET,  the 
point  must  be  on  an  inlet,  etc. 
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PART  III:  NOTATION 


10.  FORTRAN  variables  and  parameters  are  identified  below,  with 
the  corresponding  algebraic  quantities  in  the  equations  of  Thompson  and 
Bernard  (1985)  shown  in  parentheses.  The  array  dimensions  are  also 
given  in  parentheses  for  each  type  of  array. 

Main  flow  variable  field  arrays  (0: IDIM.O: JDIM)* 

VELX 

:  x  and  y  velocity  components  (u,v) 

VELY 

PRES:  difference  between  true  pressure  and  hydrostatic 
pressure  (P) 

TEMP:  temperature  (T) 

QWAL:  boundary  heat  transfer 
Additional  flow  variable  field  arrays  (0: IDIM,0: JDIM) 

UR0LD:  —  J^2(pu)n  *  -  y  (pu)n  stored  quantity  from  pre¬ 
vious  time  steps  for  second-order  time  differencing 

stored  from  previous  time  step  for  first-order  time 
differencing 

VR0LD:  as  UR0LD,  with  u  replaced  by  v 

R0LD:  as  UR0LD,  with  u  replaced  by  1  and  J/At  omitted 

ER0LD:  as  UR0LD,  with  u  replaced  by  energy  e 

UACC:  x-momentum  equation  acceleration  parameter 
VACC:  y-momemtum  equation  acceleration  parameter 
TACC:  energy  equation  acceleration  parameter 
TYPE:  point  type 

STAT:  hydrostatic  pressure  with  constant  density  RH0REF  (pH) 
WIDE:  width  (B) 


ENGY:  energy  (e) 

DIVEL:  divergence  of  velocity  (V*u) 

CMAS0:  mass  residual 

CMAS:  as  UR0LD,  with  u  replaced  by  mass  residual,  and  J 

2 

replaced  by  J  BCJD^) 

Coordinate  arrays  (0: ICIM,0: JCIM) 


cartesian  coordinates  (x,y) 


SIDE:  point  type  array  LSLIT  of  the  coordinate  code  WESC0R 


Eddy  viscosity 


CART:  artificial  viscosity 


m 


RICH:  Richardson  number  for  turbulent  eddy  viscosity  (R^ 
TURH0R:  selector  for  horizontal  turbulent  eddy  viscosity 
TURBUL:  turbulence  indicator 

C0FH0R:  specified  uniform  horizontal  turbulent  eddy  viscosity 

(Dll> 

TURVER:  selector  for  vertical  turbulent  eddy  viscosity 
C0FVER:  specified  uniform  vertical  turbulent  eddy  viscosity 
(d22) 

Run  control  parameters 

ST0RIT:  first  time  step  to  be  stored;  later,  next  step  to  be 
stored 

START:  indicates  type  of  start 
ST0INT:  interval  between  stored  time  steps 
ITERS:  maximum  number  of  iterations  allowed 
PARC0N:  accepts  partial  convergence 
STEPS:  number  of  time  steps  to  be  run 
STEP:  time  step  number  on  file 


Iteration  tolerances 


UT0L  \ 

VT0L  f 

^iteration  error  tolerances  for  u  ,  v  ,  P  ,  T 
PT0L  l 

TT0L  ) 

U N0RM:  iteration  error  norm  for  x-momentum  equation  ( | | Au | j ) 

VN0RM:  iteration  error  norm  for  y-momentum  equation  ( j | Av ( j ) 

PN0RM:  iteration  error  norm  for  pressure  equation  (||ap||) 

TN0RM:  iteration  error  norm  for  energy  equation  (|  j  AT | |) 

(IUN0RM, JUN0RM) :  location  of  maximum  x-momentum  iteration  error 
(substitute  VN0RM,  etc.,  for  other  norms) 

Indices 

D,0D:  diagonal  and  off-diagonal  elements  of  stress  and  heat 
transfer  terms 

IMAX  ) 

> field  extent 
JMAX  ) 

(I,J):  indices  of  calculation  point 
(II, JJ):  indices  for  coordinate  arrays  corresponding  to  (I,J) 
Solution  options 

ENEREQ :  allows  energy  equation  to  be  deleted 

TIMFAC (TIM0RD ,L) :  coefficients  in  time-derivative  for  order 

TIM0RD: 

TIMFAC (TIM0RD,l)fn  +  TIMFAC (TIM0RD, 2) f  n~l+  TIMFAC (TIM0RD , 3) f  ”~2 

At 


CAVFAC (CAV0RD ,L) : 


coefficients  in  extrapolation  at  concave 
corners  for  order  CAV0RD 


; 


0, 


0- 


f_  =  CAVFAC(CAV0RD,l)*(f,  +  f_)  +  CAVFAC(CAV0RD,2)*(f.  +  f . ) 


C0NVEC:  selector  for  convective  differences 


TIM0RD 

CAV0RD 

FL0BAL 

F1.0W 

STRAT 

ISTRAT 


time  derivative  order 

concave  corner  extrapolation  order 

activator  for  flow  balance 

activator  for  inflow  or  outflow 

selector  for  initial  stratification 

C-line  defining  initial  stratification 

table  of  density  defining  initial  stratification 

temperature  table  defining  initial  stratification 

elevation  table  defining  initial  stratification 


Acceleration  parameters 


UACT0L 

VACT0L 

TACT0L 


iteration  error  tolerances  for  variable  acceleration 
parameters  for  u  ,  v  ,  T  ,  respectively 


IKAX+1 


used  in  calculation  of  variable  accel¬ 
eration  parameters 


—  JMAX+1 

BACC :  boundary  temperature  acceleration  parameter 

PACC:  pressure  acceleration  parameter 

VELACC:  switch  for  variable  acceleration  parameter  for 

momentum  equations 

TEMACC:  switch  for  variable  acceleration  parameter  for  energy 

equation 


VELXAC 

VELYAC 

TEMPAC 


u  ,  v  ,  T  acceleration  parameters,  respectively, 
when  variable  parameters  are  not  calculated 


PRESAC:  pressure  acceleration  parameter  (input  value) 

B0UNAC:  boundary  temperature  acceleration  parameter 
Solution  quantities 

RH0REF:  reference  density  used  for  calculation  of  hydrostatic 
pressure  (p  ) 

GX  ) 

> gravity  vector  (g.,g~) 


FL0IN:  total  inlet  flow  rate 
FL00UT:  total  outlet  flow  rate 
DELTAT:  time  step  (At) 


DELTAT: 

VELREF 

TEMREF 

PREREF 


reference  values  of  (u,v) 


.gravity  unit  vector 


p  ,  respectively 


general  values  for  u  ,  \ 

tively,  over  entire  field 


B  ,  respec- 


(IREF.JREF): 


(XREF.YREF) : 


location  of  PREREF,  reference  pressure  for 
calculation  of  hydrostatic  pressure 


GRAC:  acceleration  of  gravity 


Output  options 


VUNITS  j 

PUNITS  >  output  units  for  (u,v)  ,  P  ,  and  T  ,  respectively 
TUN ITS  ) 

CUNITS:  output  units  for  time 
EXT0UT:  switch  for  output  of  field  extrema 

MAP0UT:  switch  for  output  of  field  maps 

MAP(5):  field  maps  desired 

RIT0UT:  switch  for  printed  output  of  field  at  selected  steps 
R1TERR:  switch  for  printed  output  of  iteration  error  norms 

RITINT:  switch  for  printed  output  of  initial  field  solution 

LABEL:  output  label 

ACCPRT:  switch  for  printed  output  concerning  variable 


ACCPRT : 


R1TEXT: 


acceleration  parameters 

switch  for  printed  output  of  outlet  variables 


Computation  site  arrays  (25) 

U  u 
V  v 
P  P 
T  T 
RH0  p 
RH0U  pu 
RH0V  pv 
RH0E  pe 
ENERGY  e 
MU  u 
KC  < 

Q  heat  source  flux 
POINT  point  type 

Sll  Ojj*  etc.  > 


SIIBAR  o^,  etc. 


Q1  qL.  etc. 
QIBAR  q ^ ,  etc. 


(25,2)  second  index  refers  to 
diagonal  (1)  or  off- 
(25,2)  diagonal  (2) 

(25,2) 

(25,2) 


DSXI 

(2) 

DSET 

(2) 

UXXI 

ux 

£ 

UYXI 

uyc 

UXET 

ux 

n 

UYET 

Uyn 

UBAR 

u 

VBAR 

V 

X 

X 

Y 

y 

XXI 

x_ 

stress  terms  in  momentum  equation,  conduction 


etc.  for  v  ...»  T  . . . ,  P 


YET  y 

n 

JCB  J  =  xv  -  x  y 

£  n  n  £ 

ALPHA  a  =  x  2  +  y  2 
n  n 


GAMMA  y  =  x^2  +  y?2 


INFCEN, JNFCEN  index  increments  for  neighboring  solution  points 
INSCEN, JNSCEN  index  increments  for  neighboring  coordinate  points 
Q  heat  flux  (q) 

Other  variables 

D11,D22  horizontal  and  vertical  eddy  viscosity 
C11.C22  horizontal  and  vertical  eddy  conductivity  (C^»C2^ 
MASDEF  mass  residual  (D) 

FL0RAT  ratio  of  inlet  to  outlet  flow  rate  before  balance 
correction 

PT0T  full  pressure  (p) 

AJAC  Jacobian  (J) 

ITACMIN , UACMAX ,  minimum,  maximum,  and  average  acceleration  parameter 
UACAVG  for  velocity  (etc.  for  V  and  T  ) 

VAN0RM  iteration  error  norm  for  x-velocity  acceleration 
parameter  (etc.  for  V  and  T  ) 

0VLEx\ 

0VELY  I 

dTFMP  \Previ°us  time  step  values  of  u  ,  v  ,  T  ,  e  ,  D 
/ respectively 
0ENGY  I 

0CMAS J 

NPT  total  number  of  field  points 
TIME  time  (t) 

TIM  time  in  output  units  (t) 

MU11,MU22  viscosity  including  turbulence  (y  +  pD^.y  +  pD£2) 
MUI122  average  of  MU11  and  MU22 


KC11,KC22  conductivity  including  turbulence 


o 


DPXI.DPET 

UC0LD.VC0LD 


(<  +  p  ^  cr<  +  P  jj  C  ) 

pressure  terms  in  momentum  equations 
velocity  values  at  previous  iterate 


DEPTH 

HEIGHT 

FULDEP 


VELH0R , VEL VER 

HORCEL, VERCEL 
RH0VER 
DMXI.DMET 

DMXIC.DMETC 

EMASDEF 


depth  to  computation  site  below  top  (d) 
height  of  computation  site  above  bottom  (h) 
top-to-bottom  distance  (d^) 

velocity  magnitude  at  computation  site  (|u|)  (etc. 
for  N  ,  S  ,  E  ,  W) 

horizontal  and  vertical  derivatives  of  velocity 

magnitude  (|u|  ,  |u|  ) 

-  x  y 

horizontal  and  vertical  cell  width  (Ax, Ay) 

vertical  derivative  of  density  (p^) 

£  and  n  convective  flux  terms  (pBfu  and  pBfv) 
(Suffix  C  indicates  coefficient  of  value  of  f  at 
computation  site.  Entire  term  is  DMXI  +  DMXIC*f  , 
etc.)  *" 

mass  residual  correction  term  (zero  for  nonconserva¬ 
tive  convective  terms,  equal  to  MASDEF  for  conserva¬ 
tive  terms) 


CENFAC  total  coefficient  of  central  value  factored  from  all 
terms 

UTEMP  intermediate  x-velocity  value  (etc.  for  v  ,  P  ,  T  ) 


(etc.  for  Y  ) 


x.x  +  y,y 
£  n  3  v  n 


xry  +  x  yF 

£7  n  T  £ 


XXIXI  x 


s 


BRTY  J(Bpt) 

DXI 

DET  (V*u) 

~  ~  n 

DXIXI  (7*u)^ 

DETET  (V*u) 

—  nn 

DX1ET  (V*u) _ 

~  ~  Cn 

DX  J(7«u) 

~  ~  x 

DY  J (V*u) 

~  ~  y 

DLAP  J2V2 (V*u) 


UXIXI  u. 


UETET  u 


UXIET  u. 


ULAP  J2V2u 
DVEL  JV  *u 


(etc.  for  V  ) 


m, 


RUX  J 


RVY  J 


RHS  BJ2V2p 


ART VC  artificial  viscosity ^CART* 
ARTCC  artificial  conductivity (cART* 


PART  IV:  OPERATION 

11.  The  code  starts  by  reading  the  namelist  BEGIN,  which  contains 
a  flag  START  which  indicates  an  initial  start,  a  continuation  of  a 
partially  converged  iteration,  or  a  continuation  of  time  steps. 

Initial  Start 


12.  For  an  original  start,  the  namelist  PARAM  is  read  next,  and 
subroutine  SETVAL  is  called  to  read  the  remainder  of  the  input.  This 
subroutine  first  reads  the  label,  the  extent  of  the  coordinate  arrays 
(IMAXC , JMAXC) ,  and  the  coordinate  point  type  array  (SIDE,  the  LSLIT  of 
WESC0R)  from  the  output  of  the  coordinate  system  code  WESC0R.  The 
extent  of  the  field  arrays  (IMAX.JMAX)  is  then  calculated,  and  the  field 
point  type  array  TYPE  is  initialized  to  0UT.  For  free  field  points, 
i.e.,  points  in  the  flow  field  and  not  on  a  boundary,  TYPE  is  then  set 
to  FIELD,  and  for  points  on  any  boundary,  TYPE  is  set  to  the  default, 
WALL,  indicating  a  no-slip  boundary.  The  coordinate  system  is  then  read 
from  the  output  of  WESC0R  into  the  coordinate  arrays  XC0R  and  YC0R. 

13.  The  boundary  types  are  designated  next  by  input  as  follows. 

For  each  boundary  segment  (outer  or  obstacle)  that  is  not  to  be  no-slip, 
namelist  INPUT  is  read  specifying  the  boundary  type  (1NLF.T,  OUTLET, 
SURFACE,  or  SLIP)  and  the  indices,  (II, Jl)  and  (12, J2),  of  the  segment 
ends.  All  points  on  that  segment  are  then  set  to  this  type  in  TYPE. 
Points  on  all  boundary  segments  not  designated  otherwise  by  input  are 
typed  no-slip  by  default,  as  noted  above,  except  that  it  is  also  pos¬ 
sible  to  change  this  default  to  slip  by  the  input.  After  all  of  the 
boundary  designations  have  been  read,  the  code  determines  if  the  seg¬ 
ments  are  bottom,  left,  right,  or  top  relative  to  the  free  field,  and 
locates  concave  and  convex  corners,  placing  the  appropriate  type  classi¬ 
fication  in  TYPE  for  each  boundary  point. 

14.  The  field  arrays  for  velocity,  pressure,  temperature,  energy, 
width,  and  wall  heat  transfer  (VELX,  VELY,  PRES,  TEMP,  ENGY,  WIDE,  and 
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QUAL)  are  initialized  to  the  default  values.  The  general  field  values 
and  the  reference  pressure  and  position  then  may  be  changed  from  the 
default  values  by  a  read  of  the  namelist  INPUT  for  each  quantity  to  be 
changed.  Units  are  read  with  the  values,  and  the  code  converts  all 
quantities  to  metric.  The  units  of  the  coordinate  system  must  be 
identified,  or  metres  will  be  used,  and  the  time  step  is  set.  Boundary 
values  on  any  boundary  segment  are  then  specified  by  a  read  of  the  name- 
list  INPUT  identifying  the  segment,  the  quantity  to  be  specified,  the 
values,  and  the  units.  The  general  field  values  remain  on  any  boundary 
segment  that  is  not  changed  in  this  manner.  Initial  density  (or  tempera¬ 
ture)  stratification  may  be  specified  either  on  a  designated  £-line  or 
with  an  input  table  of  density  (or  temperature)  versus  elevation. 

15.  Specif ication  of  temperature  on  a  boundary  segment  cause':  that 
segment  to  operate  with  a  fixed  temperature.  Specification  of  heat 
transfer  on  a  segment  causes  it  to  operate  with  that  heat  transfer.  No 
specif ication  causes  a  segment  to  operate  as  an  adiabatic  wall. 

16.  Subroutine  SETVAL  then  completes  its  setup  of  the  field  and 
boundaries  by  calculating  the  components  of  the  gravity  vector,  the 
reference  values  and  their  locations,  the  hydrostatic  pressure,  and  the 
energy  for  all  points,  and  by  initializing  the  acceleration  parameter 
arrays  to  the  default  values.  The  reference  velocity  is  the  square  root 
of  the  sum  of  the  squares  of  the  maximum  magnitudes  of  the  components  on 
the  field  and  boundary.  The  reference  temperature  and  density  are  the 
general  values  on  the  field.  The  reference  pressure  and  the  reference 
position  are  items  previously  input  with  the  general  values.  Hydro¬ 
static  pressure  is  calculated  with  the  initial  density  distribution  by 
integrating  from  the  reference  position,  where  the  full  pressure  is 
equal  to  the  reference  pressure,  with  constant  density  equal  to  the 
reference  density. 

17.  The  main  program  then  initializes  the  rest  of  the  field 
arrays,  setting  the  velocity  divergence,  DIVEL,  and  the  mass  residual 
terms,  CMAS  and  CMAS0  ,  to  zero,  and  setting  the  previous  time  flow 
variables  to  the  initial  values.  The  initial  solution  is  then  printed 
if  desired,  completing  the  initialization  of  an  original  start. 
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Continued  Iteration  Start 


18.  If  the  start  is  to  be  made  from  a  partially  converged  solu¬ 
tion,  the  read  of  namelist  BEGIN  is  followed  by  a  read  of  the  partially 
converged  solution  from  file  13.  A  read  of  the  namelist  PARAM,  so  that 
parameters  can  be  changed  if  desired,  then  completes  the  setup  for  con¬ 
tinuation  of  the  iteration. 

Continued  Time  Steps  Start 

19.  Finally,  if  the  start  is  to  begin  from  the  solution  at  some 
previous  time,  the  read  of  namelist  BEGIN  is  followed  by  a  read  of  the 
solution  at  this  time  from  file  13.  Namelist  PARAM  is  then  read,  com¬ 
pleting  the  setup  for  continuation  of  the  time  steps. 

Setup  for  Time  Step 

20.  After  any  of  these  three  possible  starting  modes,  the  conver¬ 
gence  tolerances  for  velocity,  pressure,  and  temperature  are  calculated 
from  the  reference  values.  These  tolerances  are  calculated  by  multiply¬ 
ing  the  input  tolerances  by  the  corresponding  reference  values,  except 

that  the  pressure  tolerance  is  based  on  the  maximum  of  the  reference 

2 

pressure  and  the  dynamic  pressure  pV  ,  calculated  from  the  reference 
density  and  velocity.  The  reference  values  and  tolerances  are  then 
printed,  and  the  code  is  then  ready  to  proceed  with  the  time  steps.  If 
the  step  is  the  first  time  step  from  an  initial  start,  the  code  elects 
first-order  time  difference  expressions  regardless  of  the  input  speci¬ 
fication  of  the  order.  The  input  specification  is  used  at  all  subse¬ 
quent  time  steps. 

21.  If  the  start  is  not  from  a  partially  converged  solution, 
before  proceeding  with  the  calculation  for  the  solution  at  t  +  At  ,  the 
current  values  of  velocity,  temperature,  energy,  and  the  mass  residual 
term  CMAS0  at  t  are  written  in  the  corresponding  field  arrays. 
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0VELX,  0VELY,  0TEMP,  0ENGY,  and  0CMAS.  The  inhomogeneous  terms  (UR0LD, 
VR0LD,  ER0LD,  R0LD)  and  the  mass  residual  term  CMAS  are  calculated  for 
each  FIELD  point.  If  nonconservative  convective  differences  are  used, 
the  density  is  omitted  from  UR0LD,  VR0LD,  ER0LD,  and  R0LD.  This  com¬ 
pletes  the  setup  for  the  calculation  of  the  solution  at  the  new  time 
step. 

Iteration 


Boundaries  update 

22.  The  iteration  then  proceeds.  At  each  iteration,  the  norms  are 
zeroed,  all  of  the  boundary  values  are  updated,  and  then  all  of  the 
field  values  are  updated.  In  each  case  the  updated  values  replace  the 
current  values  in  the  arrays  as  they  are  calculated,  in  accordance  with 
Gauss-Seidel  iteration.  Values  on  no-slip  or  slip  boundaries  are  up¬ 
dated  by  calling  subroutine  WALLS  from  the  appropriate  entry  point, 
B0TT0M,  T0P,  LEFT,  or  RIGHT,  relative  to  the  free  field.  (Note  that  a 
point  on  the  top  of  an  obstacle  calls  B0TT0M,  etc.)  Concave  and  convex 
corners,  whether  no-slip  or  slip,  are  treated  by  calling  the  subroutines 
CAVC0R  and  VEXC0R,  respectively.  Inlet,  outlet,  and  surface  points  are 
updated  by  calling  subroutine  INFL0,  0UTFL0,  or  SRFACE,  respectively. 

The  energy  (ENGY)  is  calculated  from  the  temperature  at  all  of  the 
boundary  points.  Consequently,  this  routine  is  called  through  the 
appropriate  entry  point,  after  each  call  to  INFL0,  0UTFL0,  or  SRFACE. 
After  the  complete  boundary  update,  the  normal  velocity  components  on 
all  outlets  or  inlets  may  be  adjusted  so  that  the  total  outflow  exactly 
equals  the  total  inflow  by  calling  subroutine  BAL0UT  or  BALIN. 

Field  update 

23.  The  updating  of  the  field  values  is  done  by  a  call  to  DIFFEQ 
at  each  FIELD  point.  After  the  field  has  been  updated,  the  error  norms 
are  checked  for  convergence  to  the  specified  tolerances.  If  convergence 
has  been  attained,  then  the  previous  time  step  values  are  retrieved,  and 
the  solution  is  written  on  the  restart  file  11.  The  solution  is  then 


printed  and  stored  on  the  solution  file  12  if  the  current  time  step  is 
designated  for  output.  The  solution  then  proceeds  to  the  next  time  step 
unless  the  prescribed  number  of  steps  has  been  executed,  in  which  case 
the  solution  stops. 

Completion 

24.  If  the  maximum  number  of  iterations  allowed  is  reached  without 
convergence,  there  are  two  possible  alternatives.  If  the  input  has  so 
specified,  the  code  will  interpret  this  as  convergence  and  proceed  with 
the  output  and  continue  to  the  next  time  step  as  above.  Otherwise,  the 
partially  converged  solution  is  written  on  the  restart  file  11.  The 
values  at  the  previous  time  step  are  retrieved  and  written  also.  In 


this  case,  the  solution  stops. 


PART  V:  SUBROUTINES 


BL0CK  Data 

25.  This  data  set  contains  the  default  values. 

Subroutine  AITKKN 

26.  This  subroutine  produces  a  new  iteration  for  a  Newton-Raphson 
iterative  solution  of  a  nonlinear  equation. 

Subroutine  CAVC0R 


27.  This  subroutine  calculates  the  temperature  on  a  concave  corner 
by  linear  or  quadratic  extrapolation  along  the  two  sides  forming  the 
corner. 


Subroutine  DIFFEQ 

28.  This  subroutine  calculates  the  values  of  the  pressure,  temper¬ 
ature,  and  velocity  components  at  each  FIELD  point.  The  routine  first 
places  the  appropriate  values  of  the  coordinates  of  the  points  surround¬ 
ing  the  point  of  calculation  in  the  arrays  X  and  Y.  The  coordinate 
derivatives  (XXI,  etc.)  are  then  calculated  where  needed  and,  from 
these,  the  Jacobian  and  a  and  y  are  calculated.  Next,  the  width, 
velocity,  temperature,  energy,  and  pressure  are  placed  where  needed. 

The  density  is  calculated  from  the  temperature. 

29.  The  viscosity,  the  turbulent  eddy  viscosity,  and  the  artifi¬ 
cial  viscosity  (if  used)  are  then  calculated  where  needed.  The  viscos¬ 
ity  is  calculated  directly  from  the  temperature.  The  turbulent  eddy 
viscosity  may  be  an  input  constant  or  may  be  calculated  from  either  of 
two  models.  The  artificial  viscosity  is  calculated  from  the  mass 
residual.  The  quantities  u  and  v  ,  the  contravariant  velocity  compo¬ 
nents,  are  the  calculated  and  the  mass  residual  is  evaluated. 


30.  The  routine  then  calls  subroutines  CH0RIN,  THERM0,  and  M0MENT 
to  calculate  the  pressure,  temperature,  and  velocity,  respectively. 
Pressure 

31.  Subroutine  CH0R1N  calculates  updated  pressure  via  the  Chorin 
scheme  (Thompson  and  Bernard  1985).  The  routine  first  evaluates  the 
mass  residual  using  first-order  one-sided  differences.  The  new  pressure 
is  set  equal  to  the  old  pressure,  minus  the  mass  residual  multiplied  by 
the  pressure  acceleration  parameter  and  the  geometric  coefficient 
derived  in  Thompson  and  Bernard  (1985).  The  change  from  the  previous 
Iteration  is  noted,  and  the  new  pressure  is  placed  in  the  pressure  field 
array  PRES.  This  routine  is  called  before  and  after  calling  M0MENT,  in 
ctuer  to  improve  convergence. 

Temperature 

32.  Subroutine  THERM0  updates  the  temperature  from  the  energy  equa¬ 
tion.  (There  is  an  input  provision  for  skipping  this  routine  if  the 
aensity  is  to  be  held  constant.)  This  routine  first  calculates  the  con¬ 
ductivity  from  the  temperature  where  needed  and  then  evaluates  the  tur¬ 
bulent  eddy  thermal  diffusivity  and  the  artificial  conductivity  (if 
used)  from  the  eddy  and  artificial  viscosities  calculated  above.  The 
heat  conduction  terms  are  then  evaluated,  with  terms  containing  the 
temperature  at  the  update  point  being  identified  as  diagonal  terms  and 
the  remaining  terms  as  off-diagonal  terms  by  the  indices  D  and  0D  , 
respectively.  Finally,  the  energy  flux  terms  are  evaluated,  using  a 
choice  of  central,  upwind,  or  ZIP  differences  (Thompson  and  Bernard 
1985),  isolating  any  terms  containing  the  energy  at  the  update  point. 
Second-order  central  differences  are  used  for  all  derivatives  except  for 
upwind  convective  terms,  which  are  first-order. 

33.  The  new  temperature  at  the  update  point  is  then  calculated 
through  a  false  position  iteration,  since  the  time  derivative  and  pos¬ 
sibly  the  energy  flux  terms  contain  the  energy  at  the  update  point, 
while  the  conduction  terms  contain  the  temperature  at  this  point.  This 
Newton-Raphson  iteration  is  done  by  using  subroutine  AITKEN.  If  vari¬ 
able  acceleration  parameters  are  being  used,  this  parameter  is  updated 
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by  calling  subroutine  TEMAC  after  the  iteration  for  the  temperature  at 
the  update  point  has  converged.  The  temperature  is  then  updated,  the 
change  from  the  previous  iteration  is  noted,  the  new  temperature  is 
placed  in  the  temperature  field  array,  TEMP,  and  the  new  energy  is  cal¬ 
culated  from  the  new  temperature. 

Velocity 

34.  Subroutine  M0MENT  calculates  the  velocity  components  from  the 
momentum  equations.  The  routine  first  evaluates  the  stress  terms,  iden¬ 
tifying  the  terms  containing  the  velocity  at  the  update  point  as  diago¬ 
nal  terms  (index  D)  and  the  remainder  as  off-diagonal  terms  (index  0D) . 
The  momentum  flux  terms  are  then  evaluated  using  a  choice  of  central, 
upwind,  or  ZIP  differences,  again  identifying  terms  containing  the 
velocity  at  the  update  point.  The  pressure  terms  are  then  evaluated. 
Second-order  central  differences  are  used  for  all  deriviarives  except 
for  pressure  gradients  and  upwind  convective  terms.  The  new  velocity  is 
calculated,  with  all  terms  involving  the  component  being  evaluated  at 
the  update  point  factored  together.  The  mass  residual  correction  term 
uD  is  omitted  if  nonconservative  convective  differences  are  used.  If 
variable  acceleration  parameters  are  to  be  used,  subroutine  VELAC  is 
then  called  to  update  this  parameter.  The  new  velocity  is  then  calcu¬ 
lated.  The  change  from  the  previous  iteration  is  noted,  and  the  new 
velocity  is  placed  in  the  velocity  field  arrays,  VELX  and  VELY. 

Subroutine  INFL0 

35.  Thi^  subroutine  calculates  the  temperature  of  an  inlet  by  set¬ 
ting  the  value  equal  to  the  adjacent  value  inside  the  field.  If  FL0W  is 
input  as  '0UTFL0W,'  the  inlet  velocity  is  extrapolated  in  the  same  way. 
The  routine  also  calculates  the  mass  flow  rate  into  the  inlet  by  summing 
the  product  of  the  appropriate  contravariant  velocity  component  and  the 
density. 


Subroutine  MAPS 


36.  This  subroutine  prints  maps  of  various  field  quantities  for 
inspection.  These  may  be  elected  by  the  input  parameters  in  the  array 
MAP: 

Point  type  map — shows  the  type  for  each  point,  i.e., 
field,  inlet,  outlet,  etc. 

]).  Velocity  map — shows  interval  contours  of  the  velocity 

magnitude  by  assigning  a  different  number  from  0  to  9  to 
each  point,  depending  on  the  ratio  of  the  velocity  from 
the  minimum  value  on  the  field  to  the  difference  between 
the  maximum  and  minimum  values  on  the  field.  Intervals  of 
0.1  from  0.0  to  1.0  are  indicated  by  the  numbers  0-9. 

£.  Temperature  map — shows  interval  contours  of  the  tempera¬ 
ture  in  the  same  manner  as  described  above  for  the 
velocity. 

cl.  Pressure  map — shows  interval  contours  of  the  difference 
between  the  pressure  and  the  hydrostatic  pressure,  as  for 
the  velocity. 

£.  Density  map — shows  interval  contours  of  the  density,  as 
for  the  velocity. 

Subroutine  0UTFL0 

37.  This  subroutine  calculates  the  temperature  of  an  outlet  by 
setting  the  value  equal  to  the  adjacent  value  inside  the  field.  If  FL0W 
is  input  as  'INFL0W,'  the  outlet  velocity  is  extrapolated  in  the  same 
way.  The  routine  also  calculates  the  mass  flow  rate  out  of  the  outlet 
in  the  manner  described  above  for  INFL0. 


Subroutine  SETVAL 


38.  This  subroutine  reads  the  input  and  sets  up  the  initial  field 
as  described  above. 


Subroutine  SRFACE 


39.  This  subroutine  calculates  the  velocity  components  and  the 
temperature  on  a  simulated  free  surface  (without  movement  of  the  sur¬ 
face)  by  setting  these  values  equal  to  the  values  at  the  adjacent  points 
inside  the  field. 


Subroutine  TEMAC 


40.  This  subroutine  calculates  the  locally  optimum  acceleration 
parameter  for  the  energy  equation.  This  parameter  depends  on  the  local 
convective  and  diffusive  terms,  and  underrelaxation  is  used  unless  dif¬ 
fusion  dominates  in  both  directions.  The  change  from  the  values  ai  the 
previous  iteration  is  noted,  and  the  new  value  is  placed  in  the  field 
array  TACC. 


Subroutine  VELAC 


41.  This  subroutine  calculates  the  locally  optimum  acceleration 
parameters  for  the  momentum  equations.  This  parameter  depends  on  the 
local  convective  and  diffusive  terms,  and  underrelaxation  is  used  unless 
diffusion  dominates  in  both  directions.  The  change  from  the  values  at 
the  previous  iteration  is  noted,  and  the  new  values  are  placed  in  the 
field  arrays,  UACC  and  VACC. 


Subroutine  VEXC0R 


42.  This  subroutine  calculates  the  temperature  and  velocity  compo¬ 
nents  at  a  convex  corner.  This  is  done  by  calling  subroutine  WALLS, 
with  the  appropriate  entry  point,  for  each  of  the  two  sides  forming  the 
corner,  and  then  averaging  the  results. 


Subroutine  WALLS 


43.  This  subroutine  calculates  the  temperature  and  pressure  on 
no-slip  and  slip  boundaries.  It  also  evaluates  the  slip  velocity  on 
slip  boundaries.  There  are  four  entry  points  for  boundaries  on  the 
B0TT0M,  T0P,  LEFT,  or  RIGHT  relative  to  the  free  field.  Note  that  the 
top  of  an  obstacle  and  the  bottom  of  the  entire  region  are  classified  as 
B0TT0M. 

44.  The  routine  first  places  the  coordinates  of  neighboring  points 
in  the  arrays  X  and  Y  and  calculates  the  coordinate  derivatives.  The 
values  of  all  quantities  needed  at  the  surrounding  points  are  placed  in 
the  appropriate  arrays,  and  the  derivatives  are  evaluated  using  second- 
order  central  differences  in  the  field  and  along  the  boundary,  except  at 
corners  where  first-order  one-sided  differences  are  used.  Such  one¬ 
sided  differences  are  also  used  for  derivatives  coming  off  the  wall. 

45.  After  the  heat  conduction  terms  are  evaluated,  either  the  new 
temp  rature  is  calculated  from  the  specified  wall  heat  transfer  or  the 
wall  heat  transfer  is  calculated  from  the  specified  wall  temperature. 

The  wall  pressure  is  evaluated  by  subroutine  CH0RIN.  If  the  boundary  is 
a  slip  boundary,  then  the  new  velocity  components  are  calculated  with 
the  vorticity  and  normal  velocity  set  to  zero. 

46.  The  new  temperature  and  velocity  components  are  calculated 
using  a  specified  uniform  acceleration  parameter,  the  change  from  the 
previous  iteration  is  noted,  and  the  new  values  are  placed  in  the  field 
arrays. 


Subroutine  XTREMA 


47.  This  subroutine  calculates  and  prints  the  extrema  of  the 
velocity,  temperature,  and  pressure. 


Subroutine  BAL0UT 


48.  This  subroutine  adjusts  the  velocities  on  all  the  outlets  to 
balance  the  inlet  and  outlet  flow  rates.  This  is  done  by  multiplying 
the  velocity  normal  to  the  outlets  by  the  ratio  of  the  inlet  flow  rate 
to  the  outlet  flow  rate.  New  velocities  are  then  calculated  (for  out¬ 
lets  only)  using  the  adjusted  normal  velocity  and  the  unchanged  tangen¬ 
tial  velocity.  Use  of  this  routine  is  an  input  option. 

Subroutine  BALIN 


49.  This  subroutine  is  the  same  as  BAL0UT,  but  is  applied  to  in¬ 
lets  instead  of  outlets.  Use  of  this  routine  is  an  input  option. 
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PART  VI:  PLOT  CODES 


50.  These  codes,  called  VECT0R  and  C0NTUR,  respectively,  plot 
velocity  vectors  at  each  point  or  selected  points  in  the  field,  and  the 
contours  of  density,  temperature,  and  pressure. 

Contour  Plots  (C0NTUR) 


51.  The  contour  plots  are  done  as  described  below.  The  code  was 

written  by  Thompson,  but  the  description  is  taken  from  Thames  (1975). 

Determination  of 
contours  in  the  4,9  plane 

52.  Let  $  =  $(4,9)  be  a  function  defined  in  the  region  D* 
possessing  continuous  second  derivatives.  Since  D*  is  closed  and 
bounded,  let  m($)  =  min  { (C » n )  |  [4»9]e  D*}  and  M($)  = 

max  {$(£,9)  |  C,^]t  D*}  .  If  i  is  a  number  such  that  m($)  <  0  < 

M($)  ,  then  we  define  the  $-contour  of  $(4,9),  C^,($)  ,  as  the  set 

CT($)  =  {[4,9]  |  [4,9]e  D*  and  $(4,9)  =  $) 

Graphically,  0^(3)  is  the  curve  created  by  the  intersection  of  the 
graph  of  $(4,9)  and  the  plane  $(4,9)  =  $  .  For  plotting  convenience, 
the  curve  is  usually  projected  onto  the  (4.9)  plane.  These  ideas  are 
illustrated  in  Figure  3. 

53.  Now  suppose  that  $(4,9)  is  known  only  in  a  discrete  fashion. 

That  is,  let  the  net  function  $  .  =  $(4.»9.)  be  known  on  the  discrete 

_  i  >  3  i  J 

set  D**  =  {[4^>9j]  I  4^  =  i-1  for  1  <  i  <  I MAX  and  =  j-1  for 

1  _<  j  _<  JMAX}  .  (The  fact  that  $  .  may  only  be  an  approximation  to 

$(4,n)  is  immaterial  to  the  current  discussion.)  If  similar  defini¬ 
tions  for  m($)  and  M($)  are  made  for  $  .  on  the  set  D**  ,  then 

f  *  J 

the  $-contour  of  $  .  ,  denoted  C  ($)  again,  can  be  defined  as 


Figure  3.  Contour  -  $  of  4>(£»n) 


The  bars  over 
be  an  element  of  D** 
of  the  net  function 
Figure  4. 


$ 


and  $  are  to  indicate  that  may  not 

and  that  <f>  is  not  necessarily  one  of  the  values 
,  The  discrete  representation  is  shown  in 

i » j 


54.  Numerically,  the  import  of  the  above  discussion  is  that 
interpolation  between  the  points  of  the  discrete  set  D**  is  required 
to  determine  C^($)  .  Consider  a  portion  of  grid  D**  as  shown  in 
Figure  5a  where  1  II  <  12  <  IMAX  and  1  <J1<J2<  JMAX  .  Each 


i-Il 


1*12  (m,n)  Cm>n  (m+l,n) 


Portion  of  grid  D** 


b.  Block  (m,n) 


Figure  5.  Sample  grid  in  set  D** 


grid  block  is  labeled  by  the  (i,j)  coordinates  of  the  lower  left-hand 
corner  of  the  block.  Block  (m,n)  is  shown  on  a  larger  scale  in 
Figure  5b. 

55.  In  order  to  improve  the  plotted  resolution,  consider  subdividing 
each  block  into  four  triangles,  as  shown.  The  value  of  <|>  at  (m  +  $, 
n  +  $)  is  taken  as  the  four-point  average 


+  $,n  +  $  ^m,n  +  ^m+l,n  +  ^m,n+l  +  ^m+l,n+l 


A  local  (C»n)  coordinate  is  affixed  to  each  grid  block  as  demonstrated  in 
Figure  5b.  In  order  to  standardize  the  interpolation  procedures,  a  local 
(£,p)  coordinate  system  is  also  placed  on  each  of  the  subtriangles  as 
illustrated  in  the  series  of  drawings  given  in  Figure  6. 


P2 


* - 


_ I 


Figure  6.  Local  coordinate  systems  for  triangles 


56.  Interpolation  is  carried  out  on  each  of  the  four  triangles  for 
each  grid  block  in  the  segment  (II  £  i  <  12,  J1  _<  j  _<  J2)  of  the  set  D** 
specified.  T.n  particular,  interpolation  is  performed  along  each  of  the 
three  sides  of  each  of  the  triangles  if  the  contour  value  <P  lies 
between  the  values  at  the  ends  of  the  sides.  Let  d'  be  the  directed 
distance  from  a  given  triangle  vertex  to  the  point  on  the  triangle  side 
where  the  contour  intersects  that  side  (denoted  by  (^)) .  This  distance 
is  illustrated  in  Figure  5b  for  triangle  and  is  defined  in  an 


analagous  manner  for  the  other  triangles.  If  is  the  value  of  $ 

*  i » j 


at  a  particular  triangle  vertex  and  4>  is  the  value  of  $  at  the 

*• 

other  end  of  the  side,  then  d'  may  be  expressed  as 


d'  -  (side  length)  (4>  -  ^ ) / (4*2  “  *1^ 

For  example,  along  side  1-2  of  triangle  d'  is  given  by 

d[-2  =  (1 . 0) ($  -  ♦Dtflfn) /(♦,,„  -  Vl.n5 

57.  Noting  that  the  sides  of  the  triangle  have  lengths  1.0,  1.0//2 
and  1.0//2,  the  contour  intersection  can  be  expressed  in  the  local 
triangle  coordinates  as 

Side  Vl 

1- 2  1-d  0 

2- 3  d/2  d/2 

3- 1  (l+d)/2  (l-d)/2 

where  d  =  (<J>  -  4)1)/($2  -  and  where  l  m  1,  2,  3,  or  4  denotes  the 

triangle  number. 

58.  Once  the  contour  intersections  have  been  determined  in  the 

local  triangle  coordinates  they  must  be  transformed  to  the  grid 

block  coordinates  (f  ,n  ).  This  is  done  in  the  conventional  fashion 

sm,n  m,n 

using  orthogonal  rotation  matrices.  If  [c  U„  ]  are  the  coordinates 

x.»  p ,  Ji.»p 

of  an  intersection  in  triangle  l  ,  then 


Note  that  up  to  three  contour  intersections  can  occur  for  each  triangle 

(i.e.,  one  on  each  side).  Finally,  the  point  [(£  )  ,  (n  )  ]  is 

m,n  Jt.p  tn.n  £,p 

transformed  to  (£,p)  coordinates  by  a  simple  linear  transformation  pro¬ 
ducing  an  elementf^,  p^]  of  the  spt  C^,($) 

Transformation  to  the  physical  plane 


59.  Since  contours  in  the  (£,p)  plane  are  of  little  interest, 

C  (*)  must  be  transformed  to  the  physical  plane.  This  is  made  possible 
through  the  use  of  the  coordinate  transformation  functions  xC^.p.) 
and  y(£^,p.)  •  Again,  interpolation  is  required  since  almost  all  ele¬ 
ments  of  C^,($)  are  not  elements  of  D**  on  which  the  discrete  func- 
ions  x(^,ii J  and  y^.p^)  are  defined.  As  illustrated  in  Figure  7, 
this  implies  a  double  linear  interpolation  must  be  performed.  If 
[E^,  n^]  denotes  an  element  of  (4> )  ,  the  first  step  is  to  locate  the 
£  and  p  values  bracketing  ^  and  "p^.  .  Denoting  these  by 

and  p.,  p.  .  as  shown  in  the  figure,  the  values  of  x.  and  y,  are 
j  j+1  k  'k 

calculated  as  follows: 


=  (n,,  -  O (X. ..  -  x4)/(p4il  -  p4)  +  x. 


for  all  k=l,2,...,N  .  Similar  expressions  are  used  to  calculate  y^  . 


Figure  7.  Interpolation  for  x  and  y 


Vectors  (VECT0R) 

60.  The  velocity  vectors  are  formed  by  simply  drawing  a  line  from 
each  point,  the  length  of  which  is  proportional  to  the  velocity  magni¬ 
tude  and  the  orientation  is  that  of  the  velocity  vector.  An  arrowhead 
is  added  to  the  tip  of  the  line. 
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APPENDIX  A:  VESSEL  INPUT  INSTRUCTIONS 


WBSSEL  Input  Initructlona 


*********************************************************************** 

* 

**********************  W  E  S  S  E  L  ********************************** 
* 

*********************************************************************** 

* 

**  INPUT  INSTRUCTIONS  *********************************************** 

* 

********************************************************************* 

* 

*  FILES  %  FOR  CONTINUATION  OF  TINE  STEPS.  OR  CONTINUATION  OK  ITERATION 

*  ON  A  PARTIALLY  CONVERGED  TINE  STEP.  THE  INITIAL  SOLUTION  IS 

*  READ  FRON  UNIT  13  UNFQRNAI TED . 

* 

*  THE  LAST  TIHE  STEF'.  OR  A  PARTIALLY  CONVERGED  TIME  STEP.  IS 

*  WRITTEN  ON  UNIT  11  UNFORMATTED  FOR  RESTART. 

* 

*  THE  TINE  STEPS  SELECTED  FOR  STORAGE  ARE  WRITTEN  UNFORMATTED 

*  ON  UNIT  12. 

* 

*  COORDINATE  SYSTEM  IS  READ  UNFORMATTED  FROM  UNIT  10. 

*  (AS  URITTEN  BY  CODE  'WESCOR'  ) 

* 

********************************************************************* 

* 

*  INPUT  X  (QUANTITIES  NOT  INDICATED  OTHERWISE  ARE  FLOAIINU  POINT.) 


NAMELIST  tBEGIN  X 

START  -  START  SELECTOR  X 

'INITIAL'  X  INITIAL  START. 

'NEXT'  X  START  FROM  CONVERGED  TIME  STEP. 

'RESTART'  X  CUNVINUE  ITERATION  ON  PARTIALLY  CONVERGED  TIME  STEP. 


NAMELIST  * P A R A M  X 

NAMELIST  ITEMS  X  (ANY  ORDER.  OMISSIONS  ALLOWED) 

(DEFAULT  VALUES  IN  PARATHENSES  AT  END) 

FLOW  -  SELECTOR  FOR  FLOW  SPECIFICATION.  ('INOUT'  ) 

'INFLOW'  X  INFLOW  SPECIFIED. 

'OUTFLOW'  X  OUTFLOW  SPECIFIED. 

'INOUT'  INFLOW  AND  OUTFLOW  SPECIFIED. 

ENEREQ  -  SELECTOR  FOR  ENERGY  EQUATION.  ('YES'  OR  'NO').  ('YES') 

('NO'  OMITS  ENERGY  EUUATIUN.  KEEPING  UNIFORM  TEMPERATURE) 

F'ARCON  -  ACCEPTOR  FOR  PARTIAL  CONVERGENCE.  ('YES'  OR  'NO').  ('NO') 

('YES'  ACCEPTS  PARTIAL  CONVERGENCE  AT  MAXI MUN  NUMBER  OF  ) 

(  ITERATIONS  AND  PROCEEDS  TO  NEXT  TIME  STEP.  ) 


STEPS  -  NUMBER  OF  TIME  STEPS. 


(INTEGER).  <1> 


STORIT  -  FIRST  STEP  TO  BE  STORED. 


(INTEGER).  (1) 


A1 


STOINT 


STRAT 


ISTRAT 


UTOL 

VTOL 

PTOL 

TTOL 


UACTOL 

VACTOL 

TACTOL 


VELXAC 

VELYAC 

PRESAC 

TEMPAC 

BOUNAC 

VEL ACC 
TEMACC 


DELTAT 

CUNITS 

VUNITS 

PUNITS 

TUNITS 


TIMORD 
C A VORD 


ITERS 


GRAX 

GRAY 

EXTOUT 

MAPOUT 


RITINT 

RITOUT 

RITEXT 


INTERVAL  BETWEEN  STORED  STEPS.  (INTEGER).  <i> 

(1  STORES  EVERY  STEP.  2  EVERY  SECOND  STEP,  ETC.) 

SELECTOR  FOR  INITIAL  STRAT I H CAT  I  UN .  ('NONE') 

'NONE'  7.  NO  STRATIFICATION. 

LINE'  7.  STRATIFICATION  ACCORDING  TO  DENS  1  I  Y  OR 
TEMPERA  fURE  INPUT  ON  XI-LINt  ISIRAT. 
'TABLE'  %  STRATIFICATION  ACCORDING  TO  INPUT  TABLE  OF 
DENSITY  OR  TEMPERATURE  VS  ELEVAI1UN. 

XI-LINE  DEFINING  INITIAL  STRATIFICATION  WHEN 
STRAT  = ' LINE '  •  (1) 

(IRRELEVANT  IF  STRAT  = ' NONE '  OR  'TABLE') 

CONVERGENCE  TOLERANCE  FOR  X-VELOCITY.  (l.E-4) 

Y-VELOCITY.  (l.E-4) 
PRESSURE.  (l.E-4) 
TEMPERATURE.  (l.E-4) 
(FRACTIONS  OF  MAXIMUM  VALUES) 


CONVERGENCE  TOLERANCE  FOR  X-MOMENTUM  ACC. PARA. 

( 1 

.E- 

4  ) 

Y-MOMENVUM  ACC. PARA. 

(  1 

.  K- 

4) 

ENERGY  ACC.F'AKA. 

(  1 

.E- 

4) 

ACCELERATION  PARAMETER  FOR  X-MOMENTUM. 

(  1  . 

) 

Y-MOMENTUM. 

(t. 

) 

CHORIN  PRESSURE  EU. 

(0. 

8) 

ENERGY. 

(1. 

> 

BOUNDARY  TEMPERATURE 

, 

( 1  . 

> 

VARIABLE  ACCELERATION  PARAMETER  ACTIVATOR  FOk  MOMENTUM.  (.TRUE.) 

ENERGY.  (.TRUE.) 

(LOGICAL).  (.TRUE.  ACTIVATES  VARIABLE  PARAMETERS.) 

TIME  STEP.  (1.  ) 

OUTPUT  UNITS  FOR  TIME.  ('SEC') 

VELOCITY.  ('MPS') 

PRESSURE.  ('PASCALS') 

TEMPERATURE »  <'K'> 

(SEE  LIST  OF  UNIT  CHOICES  BELOW.  DEFAULT  IS  METRIC  UNITS) 


ORDER  OF  TIME  DERIVATIVE.  (1) 

CORNER  EXTRAPOLATION.  (1) 

(INTEGER.  1  OR  2  ) 

MAXIMUM  NUMBER  OF  ITERATIONS.  (INTEGER).  (1) 

GRAVITY  UNIT  VECTOR  X-COMPONENT  .  (  0.) 

Y-COMPONENT.  (-1.) 


ACTIVATOR  FOR  OUTPUT  OF  FIELD  EXTREMA. 

MAPS. 


('YES'  OR  '  NO') 


(  'NO'  ) 
<  'NO'  ) 


ACTIVATOR  FOR  PRINT  OF  INITIAL  SOLUTION.  ('NO') 

TIME  STEP  SOLUTION.  ('NO') 
SOLUTION  ON  ALL  OUTLETS.  ('NO') 


RITERR 


ITERATION  NORMS 


(  'YES'  ) 


('YES'  OR  'NO') 


* 

* 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

.* 

* 

* 

* 

* 


MAP ( I  )  -  ACTIVATORS  FOR  FIELD  MAPS.  ('NONE') 

'TYPE'  X  POINT  TYPE  MAP. 

'VELOCITY'  X  VELOCITY  MAP. 

'PRES'  X  PRESSURE  MAP. 

'TEMP'  X  TEMPERATURE  MAP. 

'DENS'  X  DENSITY  MAP. 

( 'NONE'  GIVES  NO  MAPS.  ) 

CONVEC  -  SELECTOR  FOR  CONVECTIVE  DERIVATIVES.  ('UPWIND') 
CONSERVATIVE  X  'UPWIND'  OR  '  CENTRAL'  . 

NON-CONSERVATIVE  X  'NONUP'  OR  'NUNCEN'  . 

ZIP  X  'ZIP'  . 

IREF.JREF  -  INDICES  OF  ELEVATION  REFERENCE.  (1.1) 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


CART  -  COEFFICIENT  IN  ARTIFICIAL  VISCOSITY.  (0.) 

(0.0  FOR  NONE.  NOMINAL  NON-ZERO  VALUE  IS  1.0  > 

TURHOR  -  SELECTOR  FOR  HORIZONTAL  TURBULENCE.  ('NONE') 
'NONE'  X  NO  TURBULENCE. 

'CONSTANT'  X  UNIFORM  EDDY  VISCOSITY  =  COFHOR. 
'EDINGER'  X  EDINGER  MODEL. 

'KENT'  X  KENT  MODEL. 

TURVER  -  SELECTOR  FOR  VERTICAL  TURBULENCE*  ('NONE') 
'NONE'  X  NO  TURBULENCE. 

'CONSTANT'  X  UNIFORM  EDDY  VISCOSITY  =  COFVEk. 
'EDINGER'  X  EDINGER  MODEL. 

'KENT'  X  KENT  MODEL. 


*  COFHOR 

* 

* 

*  COFVER 

* 

* 

*  FLOBAL 

* 

* 

* 

* 

* 

* 

* 

* 

*  ACCPRT 

* 

* - 

* 

***  BOUNDARY  TYPE  NAMELIST  tlNPUT  X  (ANY  ORDER.  OMISSIONS  ALLOWED) 

* 

*  BEGIN  WITH  A  NAMELIST  *INPUT  WITH  ITEM= '  BOUND '  . 

* 

*  THEN  USE  A  SEPARATE  NAMELIST  FINPUT  FUR  EACH  BOUNDARY  SEGMENT. 
C 

*  CLOSE  WITH  A  NAMELIST  »INPUT  WITH  1 1  EM* ' END '  . 

* 

*  ITEM  -  BOUNDARY  TYPE  X  (DEFAULT  IS  NO-SLIP) 


-  COEFFICIENT  FOR  HORIZONTAL  TURBULENCE.  (1.) 

(IRRELEVANT  IF  TURHOR  NOT  EQUAL  'CONSTANT') 

-  COEFFICIENT  FOR  VERTICAL  TURBULENCE.  (1.) 

(IRRELEVANT  IF  TURVER  NOT  EQUAL  'CONSTANT') 

-  SELECTOR  FOR  FLOWRATE  BALANCE.  ('YES'  OR  'NO').  ('NO') 

WHEN  FLOW®' INFLOW'  THIS  CAUSES  ALL  OUTLET  NORMAL  VELOCITIES  TO 
BE  MULTIPLIED  BY  RATIO  OF  TOTAL  INLET  FLOWRATE  TO  TOTAL  OUTLET 
FLOWRATE. 

WHEN  FLOW='OUTFLOW'  THIS  CAUSES  ALL  INLET  NORMAL  VELOCITIES  TO 
BE  MULTIPLIED  RATIO  OF  TOTAL  OUTLET  FLOWRATE  TO  TOTAL  INLET 
FLOWRATE. 

WHEN  FLOW= ' INOUT '  THIS  OPTION  IS  IRRELEVANT. 

-  SELECTOR  FOR  ACCELERATION  PARAMETER  PRINT.  ('YES'  OR  'NO'). 


A3 


'  INLET ' 


'OUTLET' 


'SURFACE , 


X  INFLOW  BOUNDARY.  (VELOCITY  l  TEMPERATURE 

SPECIFIED  UR  EXTRAPOLATED) 

X  OUTFLOW  BOUNDARY.  (VELOCITY  i  TEMPERA l URE 

SPECIFIED  OR  EXTRAPOLATED) 

%  FREE  SURFACE.  (EXTRAPOLATED  VELOCITY  AND  TEMPERATURE) 
%  SLIP  BOUNDARY.  (ZERO  VURTIC1TY  I  NORMAL  VELOCITY) 


'ALL  SLIP ' %  CHANGES  BOUNDARY  DEFAULT  FROM  NO-SLIP  TO  SLIP. 
(IlfJl)  i  ( 12  >  J2)  -  INDICES  OF  ENDS  OF  BOUNDARY  SEGMENT.  (INTEGER) 


GENERAL  VALUE  NAMELIST  *1NPUT  X  (ANY  ORDER »  OMISSIONS  ALLOWED) 

BEGIN  WITH  A  NAMELIST  *INPUT  WITH  I TEM= ' GENERAL ' 

THEN  USE  A  SEPARATE  NAMELIST  *INPUT  FOR  EACH  QUANTITY. 


CLOSE  WITH  A  NAMELIST  *INPUT  WITH  I TEM=  'END ' 


ITEM  = 


' VELX '  -  X-VELOCITY  (1.) 

'VELY'  -  Y-VELOCITY  <1.) 

'PRES'  -  DIFFERENCE  FROM  REFERENCE  PRESSURE  (0.) 

'  REFPRES '  -  REFERENCE  PRESSURE  (AT  IREF.JREF)  (1  ATM) 

'TEMP'  -  TEMPERATURE  (273.15) 

'DENS'  -  DENSITY  (1000.0) 

'HEAT'  -  HEAT  TRANSFER  (0.) 

'WIDTH'  -  WIDTH  (1  .  ) 

'COORD'  -  SERVES  TO  CHANGE  COORDINATE  UNITS.  (METRIC) 
(VALUE  IS  IRRELEVANT.) 

VALUE  -  GENERAL  VALUE  OF  QUANTITY. 


UNITS  = 


(FEET/SEC) 

(MILES/HOUR) 


.  . .VELOCITY 


'KNOTS'  (NAUTICAL  MILES/HOUR) 


(ATMOSPHERES) 
(LB/SQ. IN. ) 
(LB/SQ.FT. ) 


.PRESSURE 


(CELSIUS)  .. 
(FAHRENHEIT) 


.TEMPERATURE 


IV 


CJ  * 


<  RANK INE > 


'  PCF  '  ( LB/CU  .  FT .  ) 

'  SCF  '  (SLUGS/CU.FT.  > 

'  GML '  (GRAMS /MILLILITER) 
'  GCC '  <  GRAMS/CU  •  CM • ) 


.DENSITY 


' BTUSIS '  (BTU/SQ.IN./SEC)  ....HEAT  TRANSFER 
'  BTUSFS '  (BTU/SQ.FT./SEC) 


'FEET' 

'MILES' 

'YARDS' 

'NAUTICAL  MILES' 


'FEET'  . 

'MILES' 

'YARDS' 

'NAUTICAL  MILES' 


.COORDINATES 


SECTION  VALUE  NAMELIST  ♦INPUT  %  (ANY  ORDER »  OMISSIONS  ALLOWED) 

BEGIN  WITH  A  NAMELIST  ♦INPUT  WI1H  IT EM= ' SECTI ON ' 

THEN  USE  A  SEPARATE  NAMELIST  ♦INPUT  FOR  EACH  QUANTITY 
CLOSE  WITH  A  NAMELIST  ♦INPUT  WITH  ITEM- 'END ' 

ITEM  = 


'  VELX' 


X-VELOCITY  (1.) 


'VELY'  -  Y-VELOCITY  (1.) 

'PRES'  -  DIFFERENCE  FROM  REFERENCE  PRESSURE  (1  ATM) 


'TEMP' 


TEMPERATURE  (273.15) 


DENS'  -  DENSITY  (1000.  ) 


'HEAT'  -  HEAT  TRANSFER  <0.) 


'WIDTH' 


WIDTH  <1.) 


'DENTAB'  -  INITIAL  STRATIFICATION  TABLE  X  DENSITY. 
'TEMTAB'  -  TEMPERATURE. 

'ELETAB'  -  ELEVATION. 

(IRRELEVANT  UNLESS  ' TkA I  ' f ABLE ' .  > 

(USE  'DENTAB'  OR  T.MTAB'.  NOT  BOTH.) 

VALUES  -  VALUES  OF  QUANTITY  ON  SECTION. 

<I1*J1>  i  (I2.J2)  -  INDICES  OF  CORNERS  OF  SELTIUN,  (INTEGER) 


JEL  -  NUMBER  OF  ENTRIES  IN  INITIAL  STRATIFICATION  TABLE, 
(IRRELEVANT  UNLESS  STRAT  =  '  I  ABLE ' > 


UNITS 


UNITS  OF  QUANTITY 


(SEE  LIST  ABOVE) 


APPENDIX  B:  CONTUR  INPUT  INSTRUCTIONS 


CONTUR  Input  Instructions 

jociiiiimuamttmuuiuitiituiumiiiiimitimmmmiumit 

uo-^c* 

c  c  m  t  u  r  imiwimmmmmmmmmi 

130-C* 

.  jO~C 

icO-C  CONTOUR  PLOT  CODE  -  MISSISSIPPI  STATE  UNIVERSITY  .  1932 
170-C 

:80-"C  II. S.  Ar'.rt!  ENGINEER  WATERWAYS  EXPERIMENT  STATION 

:10=C  VICKSBURG,  MISSISSIPPI 

200=C 

,’20-C 

ISO-C  InPUI  ;  i 0UAntI!IES  nCT  INDICATED  OTHERWISE  ARE  FLOATING  POINT.) 

:.o--c 

250=C«  NAMELIST  if ARAM  I 


:to=c 

nahll 

I”  ITEMS  :  An'  CRUET.  'MISSIONS  AuCWEB) 

280=0 

290=0 

300"C 

(DEFAULT  VALUES  I>  pATFNTnESES  AT  END) 

N<I 

-  number  of  r  i, <dicE;  to  be  c:an»ed.  (Integer) 

IKK 

i DEFAULT  IS  AlL 

3?0=C 

:G0-C 

mETA 

-  N.iMBEP  Of  ETA  iH'TCES  TO  BE  SCANNED.  (INTEGER' 

*<0=C 

DEFALT  IS  ALL ) 

*50=0 

ioO-C 

ISY'M 

-  SrMBOL  AnD  UNE  CONTROL.  ZERO  If  NO  SYMBOLS  ON  THE 

170-C 

CONTOUR  LINE  ARE  DESIRED.  A  VALUE  OF  3  GIVES  SYMBOLS 

330=C 

ON  EVERY  INTERSECTION,  2  OR  !  GIVES  FEWER  SYMBOLS.  A 

:oo=r 

NEGATIVE  VALUE  SUPRESSES  TtiE  CONNECTING  LINE.  SYMBOL 

•oo-c 

USED  IS  SYMBOL  »<NC0n-1,.  ( INTEGER/ 

•10=C 

: DEFAULT  IS  NO  S.MBOLS) 

oo--c 

<30  C 

NSTEPS  -  NUMBER  OF  TIME  STEFS.  (INTEGER)  (I) 

4<0=f 

*50=C 

SIZE 

•  PLOT  SIZE  In  Y-DIRECTION  IN  iMChES.  (3.) 

*oO--C 

»*o=c 

5IZRAT  -  RATIO  OF  Pl  OT  LENGTH  In  X-UIRECTICH  TO  SIZE.  (1.) 

*80=C 

»90"C 

<MIN. 

<MAX  -  CONTOUR  P' OT  LIMITS  In  X  DIRECTION. 

500=C 

(DEFAULT  IS  AlL ) 

510-  C 

520 --C 

(MIN, 

ymax  -  contour  plot  limits  in  y  direction. 

530=C 

( DEFAULT  IS  AlL  ) 

5mO--C 

C50=C 

I  OR  ID 

-  PUTS  GRID  LINES  OF  1  (ICRlDtl!  WEIGHT  AT  EACH  TIC 

5oO=C 

MARK  LOCATION.  NEGATIVE  INTEGER  PRODUCES  HALF -WEIGHT 

I'O^r 

Cm  LINES  AT  EACH  TIG  HARM  LOCATION  AnD  -JGRID  NUM¬ 

530=C 

BER  OF  FOURTH-WEIGHT  lTNCS  EVENLY  SEALED  BETWEEN) 

5=0"C 

(INTEGER)  'DEFAULT  IS  NO  GRID) 

»oo=t 

oIO-C 

DAX, 

-  HORIZONTAL  AND  VERTICAi  ACES  TIC  MARK  LENGTHS  IN  INCHES. 

o20-C 

BOX 

■  DEFAULT  IS  NO  TIC  MARKS) 

o30-C 

340=0 

OA,OB 

-  HORIZONTAL  AND  VERTICAL  AXES  MINIMUM  VALUES. 

oSO-C 

( DEFAULT  IS  NO  AXES  NUMBERS  ) 

oRO-C 

o70-C 

0A1, 

-  HORIZONTAL  AND  VERTICAL  AXES  NUMBER  INTERVAL. 

580=0 

OBI 

. DEFAULT  IS  NO  AXES  NUMBERS  ) 

590=0 

’Xl-Ctt  NAMELIST 

» INPUT  ; 

7!0--C 

*20=C 

NAMELIST  ITEMS  !  (ANY  ORDER,  OMISSIONS  ALLOWED) 

"30--C 

(A  SEPARATE  NAMELIST  1 INPUT  FOR  EACH  VARIABLE  TO 

*40=C 

750-C 

1TEM  = 

’sO=C 

"0=C 

DENSITY'  -  DENSIT; 

780=0 

790- C 

ENERGY'  -  ENERGY 

80©=f. 

SKK  PRESSURE  -  DIFFERENCE  FROM  HYDROSTATIC  PRESSURE 

820=C 

830--C  TENPERAT'  -  TEMPERATURE 

840=C 
350=0 


C  CONTOURED) 


NUM  -  NUMBER  Of  CONTOURS  FOR  THIS  VARIABLE.  ( INTEGER)  (1) 


UNITS  -- 


( DEFAULT  IS  METRIC ) 


870=C 

880=1", 


890=  C 

■prc 

-  POUNDS/ CUBIC  FOOT 

.DENS  IT: 

-•)0=c 

SCf 

-  SLDGS/CUBIf  cOCT 

;10=C 

SCC' 

-  ORAnS/CUBIC  CENTIMETER 

»20=C 

=30~C 

' FT-lBF  -  FOOT  POUNDS  . 

.ENERGY 

“♦0=C 

FT-LBi.'  -  FOOT  POUNBALS 

aSO=C 

BTU 

-  BTU 

3*0=C 

Lr^L, 

-  calories 

°'0=C 

EEC' 

-  ERGS 

?80=C 

KCAi 

-  kIi.QCAlQRIES 

3%-r 

KWH 

-  KlLOtfATT  HiOURS 

lM=C 

w* 

-  MTT  HOURS 

1010=C 

io:o=c 

'A'rt' 

-  ATrtOSPhtRES  . 

.PRESSURE 

10a0=C 

ESI 

“  P'OjK^S/SQUAal  IrtCH 

ll)40=C 

FSF" 

-  FOUNDS, SQUARE  FOOT 

1050--C 

lt)oOrC 

r  _ 

CElSIUS  . 

♦  Ttrtf ERaTuRP 

:0?0=C 

E'  - 

FAHRENHEIT 

1060=C 

R  - 

RaNKInF 

10W=C 

1100=C 

vi<t.urst  ~ 

COhIOUP  Vi'nLUt s  . 

lilO=C 

1120=8*1  IF  LESS  THAN  FOUR  VARIABLES  ARE  COnTOUKEL.  FOLLOW  Tht  CONTOUR 
ll30=C«  namELISTS  WITH  Aft  IINFUT  nahElIST  CONTAINING  ONLY  1TEH='ENU' . 

1I40=C 

1160=C 

U70=C  TIME  STEFS  ARE  REAi>  UNFORMATTED  FROM  UNIT  lit  AS  WRITTEN  BY  CODE  'NESSEl'. 

1 180=C 

1190-C  COORDINATE  system  IS  READ  UNFORMATTED  FROM  unit  IOi  AS  WRITTEN  BY  COKE  'wESCOR' 
1208=C 

ilid-cnmmuumimimmmimmmimnmmumimmu 


APPENDIX  C:  VECTOR  INPUT  INSTRUCTIONS 


VECTOR  Input  Instructions 


110=C 

:20=crowwwmwm*m  vector  mmmmmmmmmwm* 

)30=C 

loO-C  VECTOR  PLOT  CODE  -  MISSISSIPPI  STATE  UNIVERSITY  .  1982 
l70=C 

180=C  U.S.  ARMY  ENGINEER  WATERWAYS  EXPERIMENT  STATION 

l90=C  VICKSBURG;  MISSISSIPPI 

200=C 

2io-cnmnmmtnnmmtmm*mtmimnnnmnmmmtmntt 

220=C 

;30"C  SOLUTION  IS  READ  UNFORMATTED  FROM  UNIT  lb  AS  WRITTEN  BY  CODE  'WESSEL' 

240=0 

250-C  COORDINATE  SYSTEM  IS  READ  UNFORMATTED  FROM  UNIT  10;  AS  WRITTEN  BY  CODE  'WESCOR' 
260=C 

:?o=c  input  ;  (quantities  not  indicated  otherwise  are  floating  point.) 

300=C 

"iO-Ctt  NAMELIST  tPARAM  I 
320=C 

330=C  NAMELIST  ITEMS  5  (ANY  ORDER;  OMISSIONS  AllOWED) 

3+0-C  (DEFAULT  VALUES  IN  PARENTHESES  AT  END) 

350=C 

3oO=C  INSTEPS  -  TOTAl  NUMBER  OF  STEPS  TO  BE  PLOTTED.  (INTEGER)  (1) 

3SO=C  NXI  -  NUMBER  OF  XI  INDICES  TO  BE  SCANNED.  (INTEGER) 

3?0=C  (DEFAULT  IS  ALL) 

400= C 

410=C  NET A  -  NUMBER  OF  ETA  INDICES  TO  BE  SCANNED.  ( INTEGER ) 

*20=C  (DEFAULT  IS  ALL) 

+30--C 

*40=C  NSKXI  -  SKIP  INTERVAL  FOR  XI.  (INTEGER)  (1) 

450=C  (1  SCANS  EVERY  LINE.  2  EVERY  SECOND  LINE;  ETC.) 

4t»0=C 

470=C  NSKETA  -  SKIP  INTERVAL  FOR  ETA.  (INTEGER;  (1) 

480=C  (SEE  NSKXI) 

*90=0 

:00=C  SIZE  -  PLOT  SIZE  IN  Y-IURECTION  IN  INCHES.  (8.) 

510=C 

S20=C  SIZRAT  -  RATIO  OF  Pi.OT  LENGTH  IN  X-BIRECTION  TO  SIZE.  (1.) 

530=C 

S40=c  AB  -  ARROWHEAD  ANGLE  IN  DEGREES.  (20, ) 

ST;0=C 

5i>0=C  AD  -  ARROWHEAD  LENGTH  IN  INCHES.  (0.0875) 

S70=C 

530=C  S2  -  VELOCITY  SCALE  FACTOR  ,  FRACTION  OF  SIZE  PER  UNIT  VELOCITY.  (1.) 

5?0=C 

oOO-C  XMIN.XMAX  -  PLOT  LIMITS  IN  X  DIRECTION. 

5lO=C  (DEFAULT  IS  ALL) 

t>20=C 

w0=C  YMIN-YMAX  -  PLOT  LIMITS  IN  Y  DIRECTION. 

*40=C  (DEFAULT  IS  ALL) 

650=C 


APPENDIX  D:  SAMPLE  RUNSTREAMS 


WES  S  EL 

00100  /JOB 

00110  MJLQXXX<T500fCM500fP02fSTCA1 ) 

00120  USER(XXXXXXrXXXXXX) 

00130  ***  XMESSEL 

00150  FETCH <DN3fBLDrGDN-ORSBWSL«UN  =  XXXXXX»DS -Cl  ) 

00160  FETCH<DN=FT10fGDN=RSBCDOfDS=CI  > 

00170  LDR(MAP=FULLfE=1) 

00180  EXIT (U ) 

00190  COST ( LO=F  > 

00200  EXIT <U> 

00230  ST0RE(DN=FT11fGDN=RSBDT6fUN=XXXXXXfDS=CI) 

00240  EXIT( U ) 

00250  ST0RE<DN=FT12fGDN=RSBSGL6fUN  =  XXXXXXfDS  =  CI  > 

00260  EXIT ( U ) 

00270  LOGFILE<L*WESLDAY) 

00280  STORE <DN=WESLDAYfUN=XXXXXXfDS=KFfDT=C> 

00281  /EOR 

00290  E4BEGIN  START= ' INI TI AL '  ♦ 

00300  E*PARAM  STEPS=500  »  TIM0RD=2  f  PAKCON^YES'  f  IltRS--50 

00310  UT0U-1.0E-4  »  VT0L-1.0E-4  f  PTUL-l.E-4  f  TT0L“1.0E-6 

00312  STRAT* ' NONE '  f  FLOU= ' INOUT ' f  DEL1AI=5.0  * 

00321  CONVEC* ' UPWIND ' »  ENEREQ* ' YES  '  f  B0UNAC=1>0  f 

00320  PUNITS= ' ATM '  .  TUNITS='F/  .  RITEXT  = ' YES '  >  RITOUT='YES 

00330  ST0RIT=500  f  STOINT=500  t  MAPI  1  >  = 'TYPE' t  R1T INT=' YES 

00330  RITERR= ' YES '  .  PRESAC=0.B8  »  FLOBAL= ' NO  '  »  HAPUUl='YES 

00340  E4INPUT  ITEM='BOUND'  * 

00340  E4INPUT  ITEM*'SLIP'  »  11=2  .  Jl=13  f  12=16  f  J2=13  * 

00340  E* INPUT  ITEM= 'OUTLET '  ,  11  =  17  f  Jl=3  f  12  =  1/  f  J2  =  3  * 

00370  E4INPUT  ITEM=  '  END  '  * 

00380  E6INPUT  ITEM= ' GENERAL '  « 

00390  E4INPUT  ITEM= '  WIDTH '  »  VALUE=3.0  f  UNITS= ' Fthl  '  ♦ 

00400  E* INPUT  ITEM= 'VELX '  »  VALUE=0.0  f  UNITS='FPS'  * 

00400  E4INPUT  ITEM= ' VELY '  f  VALUE*0.0  f  UNITS='FPS'  * 

00410  E*INPUT  ITEM= ' TEMP '  t  VALUE=70.6  f  UNITS='F'  * 

00421  E* INPUT  I TEM= ' COORD '  f  UNITS=  'FEET  '  $ 

00430  E4INPUT  ITEM= ' END  '  * 

00440  E4INPUT  I TEM* ' SECT ION '  ♦ 

00465  E4INPUT  ITEM®' TEMP '  f  11  =  1  f  Jl=2  f  12=1  f  J2=6  f 
00465  VALUES  =  6*62.0  f  UNITS='F'  ♦ 

00481  E4INPUT  ITEM='WIDTH'  f  11=1  f  Jl=l  f  12=1  t  J2=13  f 
00481  VALUES  =  13*1.0  *  UNITS= ' FEET '  * 

00481  EtINPUT  ITEM= 'WIDTH'  *  11=2  f  Jl=l  f  12-2  ,  J2=13  f 
00481  VALUES  =  13*1.5  »  UNITS* ' FEET '  * 

00481  E*INPUT  ITEM='WIDTH'  f  11=3  .  Jl  =  l  f  12=3  f  J2  =  13  f 
00481  VALUES  =  13*2.0  f  UNITS* ' FEET'  * 

00481  E*INPUT  ITEM= ' WIDTH'  t  11=4  »  Jl  =  l  »  12=4  r  J2=13  f 
00481  VALUES  =  13*2.5  »  UNITS* ' FEET '  * 

00481  EilNPUT  I fEM= 'WIDTH '  f  11=5  »  Jl=l  »  12=5  f  J2=13  f 
00481  VALUES  =  13*3.0  »  UNITS* ' FEET '  * 

00481  E*INPUT  ITEM* ' VELX '  f  11  =  1  f  Jl*2  f  12=1  f  J2  =  7  f 
00481  VALUES  =  6*0.0446  f  UNITS='FPS'  ♦ 

00481  EtINPUT  ITEM® ' VELX '  f  11=2  f  Jl  =  2  f  12=2  f  J2= 13  f 

00481  VALUES  =  12*0.0149  f  UNITS*'FPS'  * 

00481  E*INPUT  ITEM= ' VELX '  f  11=3  r  Jl=2  f  12=3  f  J2-13  f 

00481  VALUES  =  12*0.0112  »  UNITS='FPS'  * 

00481  EBINPUT  ITEM='VELX'  f  11=4  f  Jl=2  »  12=4  f  J2=13  f 

00481  VALUES  =  12*0.0089  t  UNITS='FPS'  * 

00481  EtINPUT  ITEM* ' VELX '  t  11*5  f  Jl*2  f  12  =  5  f  J2=13  f 

00481  VALUES  =  12*0.0074  f  UNITS*'FPS'  ♦ 


00481  E*INPU  f  ITEM* ' VELX '  »  11*6  .  Jl-2  .  12  =  6  .  J2*13  » 

00481  VALUES  «  12*0.0064  .  UNITS-'FPS'  • 

00481  E*INPUT  ITEM-'VELX'  »  11*7  r  Jl  =  2  »  12=7  *  J2=13  . 

00481  VALUES  *  12*0.0056  .  UNITS-'FPS'  * 

00481  E*INPUT  I  TEN* ' VELX '  »  11  =  8  r  Jl=2  »  12=8  »  J2  =  13  » 

00481  VALUES  =  12*0.0050  >  UNITS-'FPS'  * 

00481  E*INPUT  ITEM* ' VELX '  »  11=9  .  Jl  =  2  .  12=9  »  J2=13  . 

00481  VALUES  =  12*0.0045  #  UNITS-'FPS'  * 

00481  E* INPUT  ITEM- ' VELX '  »  11-10  .  Jl  =  2  i  12*10  »  J2-13  , 
00481  VALUES  »  12*0.0041  »  UNITS-'FPS'  « 

00481  E* INPUT  ITEM* ' VELX '  ,  11=11  .  Jl=2  »  12-11  »  J2-13  . 
00481  VALUES  -  12*0.0037  t  UNITS-'FPS'  » 

00481  E*INPUT  I TEM= ' VELX '  .  11=12  .  Jl-2  f  12=12  »  J2=13  » 
00481  VALUES  -  12*0.0034  »  UNITS-'FPS'  « 

00481  E*INPUT  ITEM* '  VELX '  i  11  =  13  .  Jl-2  »  12-13  i  J2  =  13  > 
00481  VALUES  =  12*0.0032  »  UNITS-'FPS'  • 

00481  E*INPUT  ITEM= ' VELX '  .  11*14  .  Jl-2  .  12=14  ,  J2-13  » 
00481  VALUES  =  12*0.0030  t  UNITS-'FPS'  « 

00481  EtINPUT  ITEM* '  VELX '  t  11  =  15  .  Jl-2  »  12-15  .  J2-13  . 
00481  VALUES  =  12*0.0028  »  UNITS-'FPS'  * 

00481  E«INPUT  ITEM- ' VELX '  ,  11-16  .  Jl-2  »  12=16  ,  J2=13  * 
00481  VALUES  -  12*0.0026  »  UNITS*'FPS'  * 

00481  EtINPUT  ITEM- ' VELX '  *  11-17  »  Jl-3  »  12=17  ,  J2-3  » 
00481  VALUES  =  0.0297  .  UNITS-'FPS'  * 

00510  E»INPUT  ITEM*' END '  • 

00520  /EOR 
00530  /EOF 


VECTOR 


00100  /JOB 

00110  M JLQXXX<  T20»  CM300«  P02»  STCA1 ) 

00120  USER (XXXXXXf XXX XXX) 

00130  ***  XVECTOR 

00140  FETCH < DN«*BLD » GDN-OVECTOR. UN-XXXXXX f DS-CI > 

00150  FETCH ( DN-FT 10f OON-RSBCDO f OS=CI ) 

00160  FETCH (DN-FTllf  GDN-RSBS0L6 »  DS-C I ) 

00170  ACCESS(DN-CCLIB»UN=EKSAPP> 

00180  ASSIGN< DN-TAPE03 »A-FT03) 

00190  LDR (LIB-CCLIB. MAP-FULL »E«1) 

00200  EXIT (U) 

00210  COST (LO-F ) 

00220  EXIT <U> 

00250  STORE ( DN-TAPE03 * GDN-TAPE03 »UN-XXXXXX » DS=FF» DT-C > 

00260  EXIT (U) 

00270  LOGF ILE ( L-VPLTDA Y  > 

00280  STORE (DN-VPLTDAY»UN-XXXXXX*DS-FFf DT-C) 

00290  /EOR 

00300  EBPARAM  NSTEPS-1  *  S2-0.1  »  SIZE-3.0  >  SIZRAT-2.57  * 
00320  /EOF 


CONTUR 


OOIOO  /JOB 

00110  MJLQXXX<T20»CM300»P02»STCA1) 

00120  USER(XXXXXX»XXXXXX) 

00130  ***  XCONTUR 

00140  FETCH<DN=*BLD»GDN=GCQNTUR>UN=XXXXXX>DS=CI) 

00150  FErCH<BN-FT10fGl)N=RSBC00»0S-CI> 

00160  FETCH <DN=FT11,GHN=RSBS0L6»DS=CI  ) 

00170  ACCESS(DN=CCLIB>UN=£KSAPP) 

00180  ASSIGN(DN=TAPE02» A=FT02) 

00190  LDR<LIB=CCLI8,MAP*FULL»E=1> 

00200  EXI T  <  U ) 

00210  COST  <  LO  =  F ) 

00220  EXIT(U) 

00250  STORE < HN=TAPE02 » GDN=TAPE02 . UN= XXX XXX »DS=FF, UT=C > 

00260  EXIT  <  U  > 

00270  LOGFILE<L=CPLTDAY) 

00280  STORE <DN=CPLTDAY»UN=XXXXXX.DS=FF»DT=C> 

00290  /EOR 

00300  E*PARAM  NSTEPS=1  »  SIZE=3.0  .  SIZRAT=2.57  ♦ 

00360  EilNPUT  ITEH= ' TEHPERAT '  »  NUH=27  »  UNITS='F'  , 

00370  VALUES®  60.0  »  60.5  >  61.0  >  61.5  t  62.0  »  62.5  >  63.0  *  63.5  * 

00370  64.0  ,  64.5  »  65.0  >  65.5  »  66.0  >  66.5  >  67.0  >  67.5  r 

00370  68.0  »  68.5  .  69.0  .  69.5  >  70.0  »  70.5  »  71.0  i  71.5  . 

00370  72.0  .  '2.5  f  73.0  * 

00400  E4INPUT  ITEH  =  ' ENU '  * 

00420  /EOF 


I"*-- 
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APPENDIX  E:  NOTATION 


lateral  width 


mass  residual 


internal  (thermal)  energy  per  unit  mass 
arbitrary  function 
components  of  gravity  vector 

Vn  -  VC 

superscript  indicating  time-step  number 


pressure 


hydrostatic  pressure 
reference  pressure 
components  of  heat  flux  vector 


temperature 


components  of  fluid  velocity 
fluid  velocity  vector 

x-  and  y-components  of  u  ,  respectively 

fluxes  of  u  through  surfaces  of  constant  €  and  constant 
n  i  respectively 

cartesian  coordinates 

2  2 

Vyn 

Vn  +  Vn 


2  2 

xc  +  yc 


gradient  operator 


laplacian  operator 
incremental  operator 
dynamic  viscosity 
thermal  conductivity 
density 

initial  density 
reference  density 
components  of  stress  tensor 
components  of  shear  stress  tensor 
curvilinear  coordinates 


acceleration  parameter 
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